# get the command arguments from perl

# v - Input vector of values
# delta - threshold for local minima maxima
# Returns - Vector of size v, 0-nothing, 1-maxima, -1-minima
peakdet <- function (v, delta){
   x <- rep(0,length(v))
   mn = Inf
   mx = -Inf
   mnpos = NaN
   mxpos = NaN
   lookformax = 1
   for (i in 1:length(v)){
      this = v[i]
      if(this > mx){
         mx = this
         mxpos = i
      }
      if(this < mn){
         mn = this
         mnpos = i
      }
   if(lookformax==1){
    if(this < mx-delta || i==length(v)){
      x[mxpos]=1
      mn = this
      mnpos = i
      lookformax = 0
    }
   }else{
    if(this > mn+delta || i==length(v)){
      x[mnpos]=-1
      mx = this
      mxpos = i
      lookformax = 1
    }
   }
  }
  x 
}

# Input - run, runfolder, delta
args=(commandArgs(TRUE))
for(i in 1:length(args)){
    eval(parse(text=args[[i]]))
}	 
#run <- "rep1"
#runfolder <- "/home/rumen/simulation1/K=0_P=0.25_MU=1e-06_NU=1e-04_RSM=1_B=0.02_D=0.016_NN=100_T=14600_H=300_W=300"
#delta <- 0.01

dlog <-  read.csv(paste(runfolder,"/",run,".divergence.biopsy.sampling.locations.log",sep=""),header=FALSE)

peaks <- peakdet(dlog[,5],delta)
div.max.ind <- which(peaks==1)
res <- matrix(nrow=0,ncol=3)
# Include divergence of 0 for time 0, (initial "minimum")
res <- rbind(res,cbind(0,0,0))

if(length(div.max.ind)>0){
 res <- rbind(res,cbind(dlog[div.max.ind,1],dlog[div.max.ind,5],1))
}
div.min.ind <- which(peaks==-1)
if(length(div.min.ind)>0){
res <- rbind(res,cbind(dlog[div.min.ind,1],dlog[div.min.ind,5],-1))
}
if(length(which(dlog[div.min.ind,1]==max(dlog[,1])))==0&length(which(dlog[div.max.ind,1]==max(dlog[,1])))==0){
 res <- rbind(res,cbind(max(dlog[,1]),dlog[which(dlog[,1]==max(dlog[,1])),5],0))
}
colnames(res) <- c("day","divergence","max/min")
res <- res[with(as.data.frame(res),order(day)),]
write.csv(res,file=paste(runfolder,"/",run,".divergence.statistics.log",sep=""),row.names=FALSE,quote=FALSE)




